***load data

use "${path}\data\acs\acs_ready.dta", clear 

***Main Tables

cd "${path}\output\tables\"

***Table I 
*Panel C (Run NAEP and GSS Analysis before for Panel A and Panel B)

eststo: reg lifesciences evo_score if sample_acs==1, cluster(bpl) baselevels
estadd ysumm
eststo: reg lifesciences evo_score $controls  if sample_acs==1, cluster(bpl) baselevels
estadd ysumm
eststo: reg lifesciences evo_score i.bpl i.birthyr  if sample_acs==1, cluster(bpl) baselevels
estadd ysumm
eststo: reg lifesciences evo_score $controls  i.bpl i.birthyr  if sample_acs==1, cluster(bpl) baselevels
estadd ysumm
esttab using main_table_all.tex, append keep(evo_score) b(3) se(3) fragment  nomtitles nonumbers nonotes  label ///
stats(ymean ysd  r2_a N, labels ("Mean of Dep. Var." "Std. Dev. of Dep. Var." "Adj. R-squared" "Observations") fmt(2 2 3 %9.0fc))  /// 
star(* 0.1 ** 0.05 *** 0.01)  ///
refcat(evo_score "\specialcellleft{\textbf{\emph{Panel C: Outcome: Working in Life Sciences}}}", nolabel)
eststo clear

***Table II 
*Panel C (Run NAEP and GSS Analysis before for Panel A and Panel B)

eststo: reg lifesciences evo_score $controls  i.bpl i.birthyr if closestate10==1, vce(cluster bpl) baselevels
estadd ysumm
eststo: reg lifesciences evo_score $controls democrat i.bpl i.birthyr , vce(cluster bpl) baselevels
estadd ysumm
eststo: reg lifesciences evo_score $controls $statelineartrend $sq_statelineartrend i.bpl i.birthyr, vce(cluster bpl) baselevels
estadd ysumm
rename evo_score evo_score_temp
rename evo_score_dosage evo_score
eststo: reg lifesciences evo_score $controls  i.bpl i.birthyr  , cluster(bpl) baselevels
estadd ysumm
rename evo_score evo_score_dosage
rename evo_score_temp evo_score
eststo: reg lifesciences evo_score $controls i.bpl i.birthyr if onereform==1, vce(cluster bpl) baselevels
estadd ysumm
eststo: reg lifesciences evo_score $controls i.bpl i.birthyr if largestate==1, vce(cluster bpl) baselevels
estadd ysumm
eststo: reg lifesciences evo_score $controls i.bpl i.birthyr if textstate==1, vce(cluster bpl) baselevels
estadd ysumm
eststo: reg lifesciences evo_score $controls i.bpl i.birthyr if  highschoolentryyear>1994, vce(cluster bpl) baselevels
estadd ysumm
eststo: reg lifesciences evo_score $controls i.bpl i.birthyr if highschoolentryyear>1999, vce(cluster bpl) baselevels
estadd ysumm
replace evo_score = evo_score/100
probit lifesciences evo_score $controls i.bpl i.birthyr, vce(cluster bpl) baselevels
estadd ysumm
estadd scalar r2_a =round(e(r2_p), .001)
estpost margins, dydx(evo_score)
eststo ames2
eststo: reg lifesciences  $controls i.bpl i.birthyr, vce(cluster bpl) baselevels
estadd ysumm
esttab using further_robustness_all.tex, append keep(evo_score) b(3) se(3) fragment noobs nonumbers nomtitles nonotes  label ///
refcat(evo_score "\specialcellleft{\textbf{\emph{Panel C: Outcome: Working in Life Sciences}}}", nolabel) ///
star(* 0.1 ** 0.05 *** 0.01)
eststo clear

***Online Appendix Tables

cd "${path}\output\tables\"

***Table A.IV
*Panel C

centile (evo_score) if less_evolution==1, centile (20 80)

csdid lifesciences_conv if less_evolution==1 , time(highschoolentryyear) notyet gvar(first_treat) method(dripw) agg(event) rseed(2) cluster(bpl)
estat simple

csdid lifesciences_conv if less_evolution==1 & evo_score<0.95 & evo_score>0.64, time(highschoolentryyear) notyet gvar(first_treat) agg(event) method(dripw) rseed(2) cluster(state16)
estat simple

centile (evo_score) if more_evolution==1, centile (20 80)

csdid lifesciences if more_evolution==1, time(highschoolentryyear) notyet gvar(first_treat) method(dripw) agg(event) rseed(2) cluster(bpl)
estat simple

csdid lifesciences if  more_evolution==1 & evo_score<0.84 & evo_score>0.09, time(highschoolentryyear) notyet gvar(first_treat) method(dripw) agg(event) rseed(2) cluster(bpl)
estat simple

***Table A.XIV

eststo: reg lifesciences evo_score $controls  i.bpl i.birthyr if closestate10==1, vce(cluster bpl) baselevels
estadd ysumm
eststo: reg lifesciences evo_score $controls democrat i.bpl i.birthyr , vce(cluster bpl) baselevels
estadd ysumm
eststo: reg lifesciences evo_score $controls $statelineartrend $sq_statelineartrend i.bpl i.birthyr, vce(cluster bpl) baselevels
estadd ysumm
rename evo_score evo_score_temp
rename evo_score_dosage evo_score
eststo: reg lifesciences evo_score $controls  i.bpl i.birthyr  , cluster(bpl) baselevels
estadd ysumm
rename evo_score evo_score_dosage
rename evo_score_temp evo_score
eststo: reg lifesciences evo_score $controls i.bpl i.birthyr if onereform==1, vce(cluster bpl) baselevels
estadd ysumm
eststo: reg lifesciences evo_score $controls i.bpl i.birthyr if largestate==1, vce(cluster bpl) baselevels
estadd ysumm
eststo: reg lifesciences evo_score $controls i.bpl i.birthyr if textstate==1, vce(cluster bpl) baselevels
estadd ysumm
eststo: reg lifesciences evo_score $controls i.bpl i.birthyr if  highschoolentryyear>1994, vce(cluster bpl) baselevels
estadd ysumm
eststo: reg lifesciences evo_score $controls i.bpl i.birthyr if highschoolentryyear>1999, vce(cluster bpl) baselevels
estadd ysumm
replace evo_score = evo_score/100
probit lifesciences evo_score $controls i.bpl i.birthyr, vce(cluster bpl) baselevels
estadd ysumm
estadd scalar r2_a =round(e(r2_p), .001)
estpost margins, dydx(evo_score)
eststo ames2
replace evo_score = evo_score*100
esttab using acs_further_robustness.tex, replace keep(evo_score) b(3) se(3) fragment nomtitles nonotes  label ///
stats(ymean ysd r2_a N, labels ("Mean of Dep. Var." "Std. Dev. of Dep. Var." "Adj. R-squared" "Observations") fmt(2 2 3 %9.0fc))  ///
mlabels("\specialcell{Close \\ Elections}" "\specialcell{Control: \\ Governor's \\ Party}" "\specialcell{State Specific \\ Time Trends}" "\specialcell{Dosage \\ Treatment}" "\specialcell{Only One \\ Reform Event}" "\specialcell{Only Large \\ States}" "\specialcell{Only States \\ With Std. Text}""\specialcell{Sample Start: \\ 1995}" "\specialcell{Sample Start: \\ 2000}" "Probit") ///
indicate("State FEs = *bpl" "Cohort FEs = *birthyr" "Controls = female" , labels("YES" "YES")) ///
star(* 0.1 ** 0.05 *** 0.01)
eststo clear

***test for statistical significance

fvset base 5 bpl
fvset base 1990 birthyr
fvset base 2009 year

reg lifesciences evo_score $controls  i.bpl i.birthyr
estimates store base

reg lifesciences evo_score $controls  i.bpl i.birthyr if closestate10==1
estimates store close

reg lifesciences evo_score $controls democrat i.bpl i.birthyr
estimates store demo

*note: drop connecticut to avoid that variance matrix of t-test below is highly singular
reg lifesciences evo_score $controls i.bpl i.birthyr $statelineartrend $sq_statelineartrend if bpl!=9
estimates store trends

rename evo_score evo_score_temp
rename evo_score_dosage evo_score
reg lifesciences evo_score $controls  i.bpl i.birthyr
estimates store dosage
rename evo_score evo_score_dosage
rename evo_score_temp evo_score

reg lifesciences evo_score $controls  i.bpl i.birthyr if onereform==1
estimates store oneref

reg lifesciences evo_score $controls  i.bpl i.birthyr if  largestate==1
estimates store largestate

reg lifesciences evo_score $controls  i.bpl i.birthyr if textstate==1
estimates store textstate

reg lifesciences evo_score $controls  i.bpl i.birthyr if  highschoolentryyear>1994
estimates store samp1

reg lifesciences evo_score $controls  i.bpl i.birthyr if  highschoolentryyear>1999
estimates store samp2

suest base close, vce(cluster bpl) coeflegend
test _b[base_mean:evo_score] = _b[close_mean:evo_score]

suest base demo, vce(cluster bpl) coeflegend
test _b[base_mean:evo_score] = _b[demo_mean:evo_score]

suest base trends, vce(cluster bpl) coeflegend
test _b[base_mean:evo_score] = _b[trends_mean:evo_score]

suest base dosage, vce(cluster bpl) coeflegend
test _b[base_mean:evo_score] = _b[dosage_mean:evo_score]

suest base oneref, vce(cluster bpl) coeflegend
test _b[base_mean:evo_score] = _b[oneref_mean:evo_score]

suest base largestate, vce(cluster bpl) coeflegend
test _b[base_mean:evo_score] = _b[largestate_mean:evo_score]

suest base textstate, vce(cluster bpl) coeflegend
test _b[base_mean:evo_score] = _b[textstate_mean:evo_score]

suest base samp1, vce(cluster bpl) coeflegend
test _b[base_mean:evo_score] = _b[samp1_mean:evo_score]

suest base samp2, vce(cluster bpl) coeflegend
test _b[base_mean:evo_score] = _b[samp2_mean:evo_score]

quietly generate yreg = lifesciences
quietly generate yprobit = lifesciences
gen evo_score_probit = evo_score/100
gsem (yprobit  <- evo_score_probit $controls  i.bpl i.birthyr, probit  vce(cluster bpl)) (yreg <- evo_score $controls  i.bpl i.birthyr)
margins, dydx(evo_score_probit evo_score) post
margins, coeflegend
test  _b[evo_score_probit:1bn._predict]= _b[evo_score:2._predict]

***Table A.XXVI

preserve
foreach v of varlist * {
	label variable `v' `"\hspace{0.1cm} `: variable label `v''"'
	}
eststo: estpost sum evo_score lifesciences biologyscientistocc aggriscientistocc conservationscientistocc medscientistocc female white black asian native other multiple hispanic 
esttab using acs_long_descriptive_statistics.tex, replace wide cells("mean(fmt(%9.2fc) label(Mean)) sd(fmt(%9.2fc) label(\specialcell{Std. \\ Dev.})) min(fmt(%9.2fc) label(Min.)) max(fmt(%9.2fc) label(Max.))")  label booktabs nonum fragment noobs  ///
refcat(evo_score "\emph{Treatment Variable:}" lifesciences "\emph{Main Outcome - Working in life sciences:}" biologyscientistocc "\emph{Additional Outcomes - Working in subfields of life sciences:}" female "\emph{Controls:~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}", nolabel)
eststo clear
restore

***Table A.XXVIII

eststo: reg lifesciences evo_score $controls  i.bpl i.birthyr if female==1, cluster(bpl) baselevels 
estadd ysumm
eststo: reg lifesciences evo_score $controls  i.bpl i.birthyr if female==0, cluster(bpl) baselevels 
estadd ysumm
eststo: reg lifesciences evo_score $controls i.bpl i.birthyr if race==1, cluster(bpl) baselevels 
estadd ysumm
eststo: reg lifesciences evo_score $controls i.bpl i.birthyr if race==2, cluster(bpl) baselevels 
estadd ysumm
eststo: reg lifesciences evo_score $controls i.bpl i.birthyr if hispanic==1, cluster(bpl) baselevels 
estadd ysumm
esttab using acs_subgroups_individual.tex, replace keep(evo_score) b(3) se(3) fragment  nomtitles nonotes  label ///
stats(ymean ysd  r2_a N, labels ("Mean of Dep. Var." "Std. Dev. of Dep. Var." "Adj. R-squared" "Observations") fmt(2 2 5 %9.0fc)) ///
mgroups("By Gender" "By Race/Ethnicity", pattern(1 0 1 0 0) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span})) /// 
mlabels("Females" "Males" "Whites" "Blacks" "Hispanics") /// 
indicate("State FEs = *bpl" "Cohort FEs = *birthyr" "Controls = female" , labels("YES" "NO")) ///
star(* 0.1 ** 0.05 *** 0.01)
eststo clear

***test for statistical significance

reg lifesciences evo_score $controls  i.bpl i.birthyr if female==1
estimates store Female
reg lifesciences evo_score $controls  i.bpl i.birthyr if female==0
estimates store Male
reg lifesciences evo_score $controls i.bpl i.birthyr if race==1 
estimates store White
reg lifesciences evo_score $controls i.bpl i.birthyr if race==2 
estimates store Black
reg lifesciences evo_score $controls i.bpl i.birthyr if hispanic==1 
estimates store Hispanic

suest Female Male, vce(cluster bpl) coeflegend
test _b[Female_mean:evo_score] = _b[Male_mean:evo_score]

suest White Black, vce(cluster bpl) coeflegend
test _b[White_mean:evo_score] = _b[Black_mean:evo_score]

***Table A.XXXIV
foreach i of num 9/1 { 
eststo: reg lifesciences evo_score_`i' $controls  i.bpl i.birthyr, cluster(bpl) baselevels
estadd ysumm
}
esttab using acs_lifesciences_by_level.tex, replace keep(evo_score_9 evo_score_8 evo_score_7 evo_score_6 evo_score_5 evo_score_4 evo_score_3 evo_score_2 evo_score_1) b(3) se(3) fragment nomtitles nonotes  label ///
stats(ymean ysd r2_a N, labels ("Mean of Dep. Var." "Std. Dev. of Dep. Var." "Adj. R-squared" "Observations")  fmt(2 2 5 %9.0fc))  ///
mgroups("Life Sciences", pattern(1 0 0 0 0 0 0 0 0 0 0) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span})) /// 
indicate("State FEs = *bpl" "Cohort FEs = *birthyr" "Controls = female" , labels("YES" "NO"))  ///  
star(* 0.1 ** 0.05 *** 0.01)  
eststo clear

***Main Figures

cd "${path}\output\figures"

***Figure IV

reg lifesciences evo_score $controls  i.bpl i.birthyr, cluster(bpl) baselevels
estimates store Overall

reg biologyscientistocc evo_score $controls  i.bpl i.birthyr, cluster(bpl) baselevels
estimates store Biology

reg aggriscientistocc evo_score $controls  i.bpl i.birthyr, cluster(bpl) baselevels
estimates store Agriculture_and_Food

reg conservationscientistocc evo_score $controls  i.bpl i.birthyr, cluster(bpl) baselevels
estimates store Conservation_and_Forestry

reg medscientistocc evo_score $controls  i.bpl i.birthyr, cluster(bpl) baselevels
estimates store Medical_and_Other

	coefplot Biology || Agriculture_and_Food || Conservation_and_Forestry || Medical_and_Other, keep(evo_score) xline(0) bycoefs byopts(xrescale) graphregion(color(white)) bgcolor(white) bylabels("Biology (n=3,555)" "Agriculture and Food (n=1,241)" "Conservation and Forestry (n=1,047)" "Medical and Other (n=3,674)") levels(95)
      graph export "${path}\output\figures\acs_coefplot_heterogeneities.pdf", replace
	   
***Figure V
*Panel A
		  
reg lifesciences_conv pre78 pre56 pre34 zero post01 post23 post45 $controls_indiv i.bpl i.highschoolentryyear i.year if less_evolution==1, vce(cluster bpl)

testparm pre*
testparm post*

forvalue i=1(1)7 {
	local m_`i'=e(b)[1,`i']
	local v_`i'=e(V)[`i',`i']
}

matrix input mat1_twfe= (`m_1',`m_2',`m_3',`m_4',`m_5',`m_6',`m_7')
mat colnames mat1_twfe=   g_m3 g_m2 g_m1 g_0 g_1 g_2 g_3

matrix input mat3_twfe= (`v_1',`v_2',`v_3',`v_4',`v_5',`v_6',`v_7')
mat colnames mat3_twfe=  g_m3 g_m2 g_m1 g_0 g_1 g_2 g_3

event_plot mat1_twfe#mat3_twfe, stub_lag(g_#) stub_lead(g_m#) ciplottype(rcap) plottype(scatter) ///
		graph_opt(bgcolor(white) plotregion(color(white)) graphregion(color(white)) ///
		xline(0, lpattern(dash) lcolor(gs12) lwidth(thin)) legend(off) ///
		yline(0, lpattern(solid) lcolor(gs12) lwidth(thin)) legend(off) ///
		yscale(r(-0.05 0.1)) ytick(-0.05(0.05)0.1) ///
		ylabel(-0.05(0.05)0.1, labsize(small)) ///
		xscale(r(-3.5 2.5)) xtick(-3(1)2) xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small)) ///
		ytitle("Working in Life Sciences", height(6) size(small)) title("TWFE", size(small)) ///
		xtitle("Year to/from Evolution Curriculum Reform", height(6) size(small)) xsize(5) ysize(4)) ///
		lag_opt(msymbol(D) mcolor(black) msize(small)) lead_opt(msymbol(D) mcolor(black) msize(small)) ///
		lag_ci_opt(lcolor(black) lwidth(medthin)) lead_ci_opt(lcolor(black) lwidth(medthin))

graph save "event_study_TWFE.pdf", replace

csdid lifesciences_conv if less_evolution==1, time(highschoolentryyear) gvar(first_treat) agg(event) method(dripw) rseed(2) cluster(bpl)
csdid_plot
estat all

estat cevent, window(-20 -7)
local b_13 = r(table)[1,1]

estat cevent, window(-6 -5)
local b_14 = r(table)[1,1]

estat cevent, window(-4 -3)
local b_15 = r(table)[1,1]

estat cevent, window(-2 -1)
local b_16 = r(table)[1,1]

estat cevent, window(0 1)
local b_17 = r(table)[1,1]

estat cevent, window(2 3)
local b_18 = r(table)[1,1]

estat cevent, window(4 20)
local b_19 = r(table)[1,1]

matrix define mat1_cs=(`b_13',`b_14',`b_15',`b_16',`b_17',`b_18',`b_19')
mat colnames mat1_cs = T-3 T-2 T-1 T+0 T+1 T+2 T+3

estat cevent, window(-20 -7)
local v_13 = r(table)[2,1]^2

estat cevent, window(-6 -5)
local v_14 = r(table)[2,1]^2

estat cevent, window(-4 -3)
local v_15 = r(table)[2,1]^2

estat cevent, window(-2 -1)
local v_16 = r(table)[2,1]^2

estat cevent, window(0 1)
local v_17 = r(table)[2,1]^2

estat cevent, window(2 3)
local v_18 = r(table)[2,1]^2

estat cevent, window(4 20)
local v_19 = r(table)[2,1]^2

matrix define mat2_cs=(`v_13',`v_14',`v_15',`v_16',`v_17',`v_18',`v_19')
mat colnames mat2_cs = T-3 T-2 T-1 T+0 T+1 T+2 T+3

event_plot mat1_cs#mat2_cs , stub_lag(T+#) stub_lead(T-#)  ciplottype(rcap) plottype(scatter) ///
		graph_opt(bgcolor(white) plotregion(color(white)) graphregion(color(white)) ///
		xline(0, lpattern(dash) lcolor(gs12) lwidth(thin)) legend(off) ///
		yline(0, lpattern(solid) lcolor(gs12) lwidth(thin)) legend(off) ///
		yscale(r(-0.05 0.1)) ytick(-0.05(0.05)0.1) ///
		ylabel(-0.05(0.05)0.1, labsize(small)) ///
		xscale(r(-3.5 2.5)) xtick(-3(1)2) xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small)) /// 
		ytitle("Working in Life Sciences", height(6) size(small)) title("Callaway and Sant'Anna (2021)", size(small)) ///
		xtitle("Year to/from Evolution Curriculum Reform", height(6) size(small)) xsize(5) ysize(4)) ///
		lag_opt(msymbol(D) mcolor(black) msize(small)) lead_opt(msymbol(D) mcolor(black) msize(small)) ///
		lag_ci_opt(lcolor(black) lwidth(medthin)) lead_ci_opt(lcolor(black) lwidth(medthin))

graph save "event_study_CS.pdf", replace	

event_plot mat1_twfe#mat3_twfe mat1_cs#mat2_cs, stub_lag(g_# T+#) stub_lead(g_m# T-#) plottype(scatter) ciplottype(rcap) together noautolegend graph_opt( xtitle("Year to/from Evolution Curriculum Reform (two-year bins)") ytitle("Working in Life Sciences") xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small))  legend(order(1 "TWFE OLS" 3 "Callaway-Sant'Anna") rows(1) region(style(none))) xline(0, lcolor(gs8%60) lpattern(dash)) yline(0, lcolor(gs8%60) lpattern(dash)) graphregion(color(white)) bgcolor(white) ylabel(, angle(horizontal))) lag_opt1(msymbol(Oh) color(dknavy)  msize(small)) lag_ci_opt1(color(dknavy%100)) lag_opt2(msymbol(Dh) color(cranberry) msize(small)) lag_ci_opt2(color(cranberry)) lag_opt5(msymbol(Sh)  msize(large)  color(black)) lag_ci_opt5(color(black)) perturb(-0.05(0.07)0.01) alpha(0.05)
graph export "${path}\output\figures\acs-less_twfecs.pdf", replace

***Figure V
*Panel B

reg lifesciences_conv pre78 pre56 pre34 zero post01 post23 post45 $controls_indiv i.bpl i.highschoolentryyear i.year if more_evolution==1, vce(cluster bpl)

testparm pre*
testparm post*

forvalue i=1(1)7 {
	local m_`i'=e(b)[1,`i']
	local v_`i'=e(V)[`i',`i']
}

matrix input mat1_twfe= (`m_1',`m_2',`m_3',`m_4',`m_5',`m_6',`m_7')
mat colnames mat1_twfe=   g_m3 g_m2 g_m1 g_0 g_1 g_2 g_3

matrix input mat3_twfe= (`v_1',`v_2',`v_3',`v_4',`v_5',`v_6',`v_7')
mat colnames mat3_twfe=  g_m3 g_m2 g_m1 g_0 g_1 g_2 g_3

event_plot mat1_twfe#mat3_twfe, stub_lag(g_#) stub_lead(g_m#) ciplottype(rcap) plottype(scatter) ///
		graph_opt(bgcolor(white) plotregion(color(white)) graphregion(color(white)) ///
		xline(0, lpattern(dash) lcolor(gs12) lwidth(thin)) legend(off) ///
		yline(0, lpattern(solid) lcolor(gs12) lwidth(thin)) legend(off) ///
		yscale(r(-0.05 0.1)) ytick(-0.05(0.05)0.1) ///
		ylabel(-0.05(0.05)0.1, labsize(small)) ///
		xscale(r(-3.5 2.5)) xtick(-3(1)2) xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small)) ///
		ytitle("Working in Life Sciences", height(6) size(small)) title("TWFE", size(small)) ///
		xtitle("Year to/from Evolution Curriculum Reform", height(6) size(small)) xsize(5) ysize(4)) ///
		lag_opt(msymbol(D) mcolor(black) msize(small)) lead_opt(msymbol(D) mcolor(black) msize(small)) ///
		lag_ci_opt(lcolor(black) lwidth(medthin)) lead_ci_opt(lcolor(black) lwidth(medthin))

graph save "event_study_TWFE.pdf", replace

csdid lifesciences_conv if more_evolution==1, time(highschoolentryyear) gvar(first_treat) agg(event) method(dripw) rseed(2) cluster(bpl)
csdid_plot
estat all

estat cevent, window(-20 -7)
local b_13 = r(table)[1,1]

estat cevent, window(-6 -5)
local b_14 = r(table)[1,1]

estat cevent, window(-4 -3)
local b_15 = r(table)[1,1]

estat cevent, window(-2 -1)
local b_16 = r(table)[1,1]

estat cevent, window(0 1)
local b_17 = r(table)[1,1]

estat cevent, window(2 3)
local b_18 = r(table)[1,1]

estat cevent, window(4 20)
local b_19 = r(table)[1,1]


matrix define mat1_cs=(`b_13',`b_14',`b_15',`b_16',`b_17',`b_18',`b_19')
mat colnames mat1_cs = T-3 T-2 T-1 T+0 T+1 T+2 T+3

estat cevent, window(-20 -7)
local v_13 = r(table)[2,1]^2

estat cevent, window(-6 -5)
local v_14 = r(table)[2,1]^2

estat cevent, window(-4 -3)
local v_15 = r(table)[2,1]^2

estat cevent, window(-2 -1)
local v_16 = r(table)[2,1]^2

estat cevent, window(0 1)
local v_17 = r(table)[2,1]^2

estat cevent, window(2 3)
local v_18 = r(table)[2,1]^2

estat cevent, window(4 20)
local v_19 = r(table)[2,1]^2

matrix define mat2_cs=(`v_13',`v_14',`v_15',`v_16',`v_17',`v_18',`v_19')
mat colnames mat2_cs = T-3 T-2 T-1 T+0 T+1 T+2 T+3


event_plot mat1_cs#mat2_cs , stub_lag(T+#) stub_lead(T-#)  ciplottype(rcap) plottype(scatter) ///
		graph_opt(bgcolor(white) plotregion(color(white)) graphregion(color(white)) ///
		xline(0, lpattern(dash) lcolor(gs12) lwidth(thin)) legend(off) ///
		yline(0, lpattern(solid) lcolor(gs12) lwidth(thin)) legend(off) ///
		yscale(r(-0.05 0.1)) ytick(-0.05(0.05)0.1) ///
		ylabel(-0.05(0.05)0.1, labsize(small)) ///
		xscale(r(-3.5 2.5)) xtick(-3(1)2) xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small)) /// 
		ytitle("Working in Life Sciences", height(6) size(small)) title("Callaway and Sant'Anna (2021)", size(small)) ///
		xtitle("Year to/from Evolution Curriculum Reform", height(6) size(small)) xsize(5) ysize(4)) ///
		lag_opt(msymbol(D) mcolor(black) msize(small)) lead_opt(msymbol(D) mcolor(black) msize(small)) ///
		lag_ci_opt(lcolor(black) lwidth(medthin)) lead_ci_opt(lcolor(black) lwidth(medthin))

graph save "event_study_CS.pdf", replace	

event_plot mat1_twfe#mat3_twfe mat1_cs#mat2_cs, stub_lag(g_# T+#) stub_lead(g_m# T-#) plottype(scatter) ciplottype(rcap) together noautolegend graph_opt( xtitle("Year to/from Evolution Curriculum Reform (two-year bins)") ytitle("Working in Life Sciences") xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small))  legend(order(1 "TWFE OLS" 3 "Callaway-Sant'Anna") rows(1) region(style(none))) xline(0, lcolor(gs8%60) lpattern(dash)) yline(0, lcolor(gs8%60) lpattern(dash)) graphregion(color(white)) bgcolor(white) ylabel(, angle(horizontal))) lag_opt1(msymbol(Oh) color(dknavy)  msize(small)) lag_ci_opt1(color(dknavy%100)) lag_opt2(msymbol(Dh) color(cranberry) msize(small)) lag_ci_opt2(color(cranberry)) lag_opt5(msymbol(Sh)  msize(large)  color(black)) lag_ci_opt5(color(black)) perturb(-0.05(0.07)0.01) alpha(0.05)
graph export "${path}\output\figures\acs-more_twfecs.pdf", replace

***Online Appendix Figures

cd "${path}\output\figures"

***Figure A.VI 

preserve
clear
cd "${path}\data\acs\"

set obs 38
set seed 2
  
  foreach num of numlist 1/1000 {
  gen reformyear`num'=.
  replace reformyear`num' = 2005 in 1
  replace reformyear`num' = 2006 in 2
  replace reformyear`num' = 2005 in 3
  replace reformyear`num' = 2009 in 4
  replace reformyear`num' = 2004 in 5
  replace reformyear`num' = 2006 in 6
  replace reformyear`num' = 2006 in 7
  replace reformyear`num' = 2008 in 8
  replace reformyear`num' = 2004 in 9
  replace reformyear`num' = 2005 in 10
  replace reformyear`num' = 2004 in 11
  replace reformyear`num' = 2006 in 12
  replace reformyear`num' = 2007 in 13
  replace reformyear`num' = 2007 in 14
  replace reformyear`num' = 2005 in 15
  replace reformyear`num' = 2007 in 16
  replace reformyear`num' = 2002 in 17
  replace reformyear`num' = 2006 in 18
  replace reformyear`num' = 2000 in 19
  replace reformyear`num' = 2009 in 20
  replace reformyear`num' = 2008 in 21
  replace reformyear`num' = 2008 in 22
  replace reformyear`num' = 2006 in 23
  replace reformyear`num' = 2004 in 24
  replace reformyear`num' = 2006 in 25
  replace reformyear`num' = 2003 in 26
  replace reformyear`num' = 2004 in 27
  replace reformyear`num' = 2006 in 28
  replace reformyear`num' = 2006 in 29  
  replace reformyear`num' = 2002 in 30
  replace reformyear`num' = 2006 in 31
  replace reformyear`num' = 2005 in 32
  replace reformyear`num' = 2005 in 33
  replace reformyear`num' = 2007 in 34
  replace reformyear`num' = 2009 in 35
  replace reformyear`num' = 2003 in 36
  replace reformyear`num' = 2008 in 37
  replace reformyear`num' = 2003 in 38
 
 gen rand`num' = runiform()
 sort rand`num'
 drop rand`num'
 }

 gen state=.
 replace state= 1 in 1
 replace state= 2 in 2
 replace state= 5 in 3
 replace state= 8 in 4
 replace state= 9 in 5
 replace state= 10 in 6
 replace state= 11 in 7
 replace state= 12 in 8
 replace state= 13 in 9
 replace state= 15 in 10 
 replace state= 17 in 11
 replace state= 18 in 12
 replace state= 19 in 13
 replace state= 20 in 14
 replace state= 22 in 15
 replace state= 23 in 16
 replace state= 24 in 17
 replace state= 25 in 18
 replace state= 26 in 19
 replace state= 27 in 20
 replace state= 28 in 21
 replace state= 29 in 22
 replace state= 30 in 23
 replace state= 32 in 24
 replace state= 33 in 25
 replace state= 35 in 26
 replace state= 37 in 27
 replace state= 38 in 28
 replace state= 39 in 29
 replace state= 42 in 30
 replace state= 44 in 31 
 replace state= 45 in 32
 replace state= 46 in 33
 replace state= 47 in 34
 replace state= 48 in 35
 replace state= 51 in 36
 replace state= 54 in 37
 replace state= 56 in 38
 
order state
save "permutation.dta", replace 
restore

merge m:1 state using "permutation.dta" 

tempname myresults 
postfile `myresults' threshold intercept gradient se using myresults.dta, replace 
 foreach i of num 1/1000 {
 
gen evolution_score_perm_`i'=.
replace evolution_score_perm_`i'=9 if state==1 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=21 if state==1 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=48 if state==2 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=59 if state==2 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=82 if state==4 
replace evolution_score_perm_`i'=55 if state==5 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=66 if state==5 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=100 if state==6
replace evolution_score_perm_`i'=86 if state==8 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=82 if state==8 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=100 if state==9 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=59 if state==9 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=91 if state==10 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=80 if state==10 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=80 if state==11 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=96 if state==11 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=16 if state==12 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=91 if state==12 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=7 if state==13 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=66 if state==13 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=91 if state==15 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=75 if state==15 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=82 if state==16 
replace evolution_score_perm_`i'=45 if state==17 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=82 if state==17 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=100 if state==18 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=96 if state==18 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=77 if state==19 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=0 if state==20 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=96 if state==20 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=55 if state==21 
replace evolution_score_perm_`i'=64 if state==22 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=27 if state==22 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=30 if state==23 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=68 if state==23 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=77 if state==24 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=73 if state==24 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=82 if state==25 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=84 if state==25 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=84 if state==26 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=80 if state==26 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=86 if state==27 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=89 if state==27 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=5 if state==28 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=87 if state==28 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=82 if state==29 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=78 if state==29 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=82 if state==30 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=75 if state==30 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=66 if state==31
replace evolution_score_perm_`i'=70 if state==32 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=77 if state==32 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=23 if state==33 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=91 if state==33 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=100 if state==34 
replace evolution_score_perm_`i'=73 if state==35 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=91 if state==35 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=68 if state==36
replace evolution_score_perm_`i'=100 if state==37 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=82 if state==37 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=9 if state==38 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=64 if state==38 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=30 if state==39 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=86 if state==39 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=25 if state==40 
replace evolution_score_perm_`i'=82 if state==41 
replace evolution_score_perm_`i'=91 if state==42 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=96 if state==42 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=100 if state==44 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=82 if state==44 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=95 if state==45 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=91 if state==45 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=82 if state==46 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=77 if state==46 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=2 if state==47 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=55 if state==47 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=64 if state==48 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=46 if state==48 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=82 if state==49 
replace evolution_score_perm_`i'=86 if state==50 
replace evolution_score_perm_`i'=50 if state==51 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=68 if state==51 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=86 if state==53 
replace evolution_score_perm_`i'=2 if state==54 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=46 if state==54 & highschoolentryyear>=reformyear`i'
replace evolution_score_perm_`i'=55 if state==55
replace evolution_score_perm_`i'=36 if state==56 & highschoolentryyear<reformyear`i'
replace evolution_score_perm_`i'=61 if state==56 & highschoolentryyear>=reformyear`i'
gen evo_score_perm_`i' = evolution_score_perm_`i'/100
label variable evo_score_perm_`i' "Evolution Score"

reg lifesciences evo_score_perm_`i' $controls  i.bpl i.birthyr  , cluster(bpl) baselevels
post `myresults' (`i') (`=_b[_cons]') (`=_b[evo_score_perm_`i']') (`=_se[evo_score_perm_`i']') 
}
postclose `myresults' 

preserve
use "${path}\data\gss\myresults.dta", clear
twoway kdensity gradient, lcolor(black) xline(0.0354, lcolor(red)) xline(0.0344, lcolor(black) lpattern(dash)) graphregion(fcolor(white)) bgcolor(white) xtitle("Coefficients from placebo regressions") ytitle("Density")
graph export "${path}\output\figures\permutations_acs.pdf", replace
sum gradient, d
restore

***Figure A.IX
*Panel A

reghdfe lifesciences_conv pre78 pre56 pre34 zero post01 post23 post45 $controls_indiv if less_evolution==1, absorb(bpl highschoolentryyear year) cluster(bpl) noconstant
matrix l_vec = 0 \ 1 \ 0
local plotopts xtitle(Mbar) ytitle(95% Robust CI)
honestdid, l_vec(l_vec) pre(1/3) post(4/6) mvec(0(0.001)0.01) omit alpha(0.05) delta(sd) coefplot `plotopts' yscale(r(-0.05 0.15)) ytick(-0.05(0.05)0.15) ylabel(-0.05(0.05)0.15)
graph export "${path}\output\figures\acs-less-rr.pdf", replace

***Figure A.IX
*Panel B

reghdfe lifesciences_conv pre78 pre56 pre34 zero post01 post23 post45 $controls_indiv if more_evolution==1, absorb(bpl highschoolentryyear year) cluster(bpl) noconstant
matrix l_vec = 0 \ 1 \ 0
local plotopts xtitle(Mbar) ytitle(95% Robust CI)
honestdid, l_vec(l_vec) pre(1/3) post(4/6) mvec(0(0.001)0.01) omit alpha(0.05)  delta(sd) coefplot `plotopts' yscale(r(-0.05 0.05)) ytick(-0.05(0.05)0.05) ylabel(-0.05(0.05)0.05)
graph export "${path}\output\figures\acs-more-rr.pdf", replace

***Figure A.XII
*Panel A

reg lifesciences_conv pre78 pre56 pre34 zero post01 post23 post45 $controls_indiv i.bpl i.highschoolentryyear i.year if less_evolution==1 | bpl==31 | bpl==41 | bpl==49 | bpl==55, vce(cluster bpl)

testparm pre*
testparm post*

forvalue i=1(1)7 {
	local m_`i'=e(b)[1,`i']
	local v_`i'=e(V)[`i',`i']
}

matrix input mat1_twfe= (`m_1',`m_2',`m_3',`m_4',`m_5',`m_6',`m_7')
mat colnames mat1_twfe=   g_m3 g_m2 g_m1 g_0 g_1 g_2 g_3

matrix input mat3_twfe= (`v_1',`v_2',`v_3',`v_4',`v_5',`v_6',`v_7')
mat colnames mat3_twfe=  g_m3 g_m2 g_m1 g_0 g_1 g_2 g_3

event_plot mat1_twfe#mat3_twfe, stub_lag(g_#) stub_lead(g_m#) ciplottype(rcap) plottype(scatter) ///
		graph_opt(bgcolor(white) plotregion(color(white)) graphregion(color(white)) ///
		xline(0, lpattern(dash) lcolor(gs12) lwidth(thin)) legend(off) ///
		yline(0, lpattern(solid) lcolor(gs12) lwidth(thin)) legend(off) ///
		yscale(r(-0.05 0.1)) ytick(-0.05(0.05)0.1) ///
		ylabel(-0.05(0.05)0.1, labsize(small)) ///
		xscale(r(-3.5 2.5)) xtick(-3(1)2) xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small)) ///
		ytitle("Working in Life Sciences", height(6) size(small)) title("TWFE", size(small)) ///
		xtitle("Year to/from Evolution Curriculum Reform", height(6) size(small)) xsize(5) ysize(4)) ///
		lag_opt(msymbol(D) mcolor(black) msize(small)) lead_opt(msymbol(D) mcolor(black) msize(small)) ///
		lag_ci_opt(lcolor(black) lwidth(medthin)) lead_ci_opt(lcolor(black) lwidth(medthin))

graph save "event_study_TWFE.pdf", replace

csdid lifesciences_conv if less_evolution==1 | bpl==31 | bpl==41 | bpl==49 | bpl==55 , time(highschoolentryyear) notyet gvar(first_treat) agg(event) method(dripw) rseed(2) cluster(bpl)
csdid_plot
estat all

estat cevent, window(-20 -7)
local b_13 = r(table)[1,1]

estat cevent, window(-6 -5)
local b_14 = r(table)[1,1]

estat cevent, window(-4 -3)
local b_15 = r(table)[1,1]

estat cevent, window(-2 -1)
local b_16 = r(table)[1,1]

estat cevent, window(0 1)
local b_17 = r(table)[1,1]

estat cevent, window(2 3)
local b_18 = r(table)[1,1]

estat cevent, window(4 20)
local b_19 = r(table)[1,1]

matrix define mat1_cs=(`b_13',`b_14',`b_15',`b_16',`b_17',`b_18',`b_19')
mat colnames mat1_cs = T-3 T-2 T-1 T+0 T+1 T+2 T+3

estat cevent, window(-20 -7)
local v_13 = r(table)[2,1]^2

estat cevent, window(-6 -5)
local v_14 = r(table)[2,1]^2

estat cevent, window(-4 -3)
local v_15 = r(table)[2,1]^2

estat cevent, window(-2 -1)
local v_16 = r(table)[2,1]^2

estat cevent, window(0 1)
local v_17 = r(table)[2,1]^2

estat cevent, window(2 3)
local v_18 = r(table)[2,1]^2

estat cevent, window(4 20)
local v_19 = r(table)[2,1]^2

matrix define mat2_cs=(`v_13',`v_14',`v_15',`v_16',`v_17',`v_18',`v_19')
mat colnames mat2_cs = T-3 T-2 T-1 T+0 T+1 T+2 T+3

event_plot mat1_cs#mat2_cs , stub_lag(T+#) stub_lead(T-#)  ciplottype(rcap) plottype(scatter) ///
		graph_opt(bgcolor(white) plotregion(color(white)) graphregion(color(white)) ///
		xline(0, lpattern(dash) lcolor(gs12) lwidth(thin)) legend(off) ///
		yline(0, lpattern(solid) lcolor(gs12) lwidth(thin)) legend(off) ///
		yscale(r(-0.05 0.1)) ytick(-0.05(0.05)0.1) ///
		ylabel(-0.05(0.05)0.1, labsize(small)) ///
		xscale(r(-3.5 2.5)) xtick(-3(1)2) xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small)) /// 
		ytitle("Working in Life Sciences", height(6) size(small)) title("Callaway and Sant'Anna (2021)", size(small)) ///
		xtitle("Year to/from Evolution Curriculum Reform", height(6) size(small)) xsize(5) ysize(4)) ///
		lag_opt(msymbol(D) mcolor(black) msize(small)) lead_opt(msymbol(D) mcolor(black) msize(small)) ///
		lag_ci_opt(lcolor(black) lwidth(medthin)) lead_ci_opt(lcolor(black) lwidth(medthin))

graph save "event_study_CS.pdf", replace	

event_plot mat1_twfe#mat3_twfe mat1_cs#mat2_cs, stub_lag(g_# T+#) stub_lead(g_m# T-#) plottype(scatter) ciplottype(rcap) together noautolegend graph_opt( xtitle("Year to/from Evolution Curriculum Reform (two-year bins)") ytitle("Working in Life Sciences") xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small))  legend(order(1 "TWFE OLS" 3 "Callaway-Sant'Anna") rows(1) region(style(none))) xline(0, lcolor(gs8%60) lpattern(dash)) yline(0, lcolor(gs8%60) lpattern(dash)) graphregion(color(white)) bgcolor(white) ylabel(, angle(horizontal))) lag_opt1(msymbol(Oh) color(dknavy)  msize(small)) lag_ci_opt1(color(dknavy%100)) lag_opt2(msymbol(Dh) color(cranberry) msize(small)) lag_ci_opt2(color(cranberry)) lag_opt5(msymbol(Sh)  msize(large)  color(black)) lag_ci_opt5(color(black)) perturb(-0.05(0.07)0.01) alpha(0.05)
graph export "${path}\output\figures\acs-less_twfecs_nevert.pdf", replace

***Figure A.XII
*Panel B

reg lifesciences_conv pre78 pre56 pre34 zero post01 post23 post45 $controls_indiv i.bpl i.highschoolentryyear i.year if more_evolution==1 | bpl==31 | bpl==41 | bpl==49 | bpl==55, vce(cluster bpl)

testparm pre*
testparm post*

forvalue i=1(1)7 {
	local m_`i'=e(b)[1,`i']
	local v_`i'=e(V)[`i',`i']
}

matrix input mat1_twfe= (`m_1',`m_2',`m_3',`m_4',`m_5',`m_6',`m_7')
mat colnames mat1_twfe=   g_m3 g_m2 g_m1 g_0 g_1 g_2 g_3

matrix input mat3_twfe= (`v_1',`v_2',`v_3',`v_4',`v_5',`v_6',`v_7')
mat colnames mat3_twfe=  g_m3 g_m2 g_m1 g_0 g_1 g_2 g_3

event_plot mat1_twfe#mat3_twfe, stub_lag(g_#) stub_lead(g_m#) ciplottype(rcap) plottype(scatter) ///
		graph_opt(bgcolor(white) plotregion(color(white)) graphregion(color(white)) ///
		xline(0, lpattern(dash) lcolor(gs12) lwidth(thin)) legend(off) ///
		yline(0, lpattern(solid) lcolor(gs12) lwidth(thin)) legend(off) ///
		yscale(r(-0.05 0.1)) ytick(-0.05(0.05)0.1) ///
		ylabel(-0.05(0.05)0.1, labsize(small)) ///
		xscale(r(-3.5 2.5)) xtick(-3(1)2) xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small)) ///
		ytitle("Working in Life Sciences", height(6) size(small)) title("TWFE", size(small)) ///
		xtitle("Year to/from Evolution Curriculum Reform", height(6) size(small)) xsize(5) ysize(4)) ///
		lag_opt(msymbol(D) mcolor(black) msize(small)) lead_opt(msymbol(D) mcolor(black) msize(small)) ///
		lag_ci_opt(lcolor(black) lwidth(medthin)) lead_ci_opt(lcolor(black) lwidth(medthin))

graph save "event_study_TWFE.pdf", replace

csdid lifesciences_conv if more_evolution==1 | bpl==31 | bpl==41 | bpl==49 | bpl==55 , time(highschoolentryyear) notyet gvar(first_treat) agg(event) method(dripw) rseed(2) cluster(bpl)
csdid_plot
estat all

estat cevent, window(-20 -7)
local b_13 = r(table)[1,1]

estat cevent, window(-6 -5)
local b_14 = r(table)[1,1]

estat cevent, window(-4 -3)
local b_15 = r(table)[1,1]

estat cevent, window(-2 -1)
local b_16 = r(table)[1,1]

estat cevent, window(0 1)
local b_17 = r(table)[1,1]

estat cevent, window(2 3)
local b_18 = r(table)[1,1]

estat cevent, window(4 20)
local b_19 = r(table)[1,1]

matrix define mat1_cs=(`b_13',`b_14',`b_15',`b_16',`b_17',`b_18',`b_19')
mat colnames mat1_cs = T-3 T-2 T-1 T+0 T+1 T+2 T+3

estat cevent, window(-20 -7)
local v_13 = r(table)[2,1]^2

estat cevent, window(-6 -5)
local v_14 = r(table)[2,1]^2

estat cevent, window(-4 -3)
local v_15 = r(table)[2,1]^2

estat cevent, window(-2 -1)
local v_16 = r(table)[2,1]^2

estat cevent, window(0 1)
local v_17 = r(table)[2,1]^2

estat cevent, window(2 3)
local v_18 = r(table)[2,1]^2

estat cevent, window(4 20)
local v_19 = r(table)[2,1]^2

matrix define mat2_cs=(`v_13',`v_14',`v_15',`v_16',`v_17',`v_18',`v_19')
mat colnames mat2_cs = T-3 T-2 T-1 T+0 T+1 T+2 T+3

event_plot mat1_cs#mat2_cs , stub_lag(T+#) stub_lead(T-#)  ciplottype(rcap) plottype(scatter) ///
		graph_opt(bgcolor(white) plotregion(color(white)) graphregion(color(white)) ///
		xline(0, lpattern(dash) lcolor(gs12) lwidth(thin)) legend(off) ///
		yline(0, lpattern(solid) lcolor(gs12) lwidth(thin)) legend(off) ///
		yscale(r(-0.05 0.1)) ytick(-0.05(0.05)0.1) ///
		ylabel(-0.05(0.05)0.1, labsize(small)) ///
		xscale(r(-3.5 2.5)) xtick(-3(1)2) xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small)) /// 
		ytitle("Working in Life Sciences", height(6) size(small)) title("Callaway and Sant'Anna (2021)", size(small)) ///
		xtitle("Year to/from Evolution Curriculum Reform", height(6) size(small)) xsize(5) ysize(4)) ///
		lag_opt(msymbol(D) mcolor(black) msize(small)) lead_opt(msymbol(D) mcolor(black) msize(small)) ///
		lag_ci_opt(lcolor(black) lwidth(medthin)) lead_ci_opt(lcolor(black) lwidth(medthin))

graph save "event_study_CS.pdf", replace	

event_plot mat1_twfe#mat3_twfe mat1_cs#mat2_cs, stub_lag(g_# T+#) stub_lead(g_m# T-#) plottype(scatter) ciplottype(rcap) together noautolegend graph_opt( xtitle("Year to/from Evolution Curriculum Reform (two-year bins)") ytitle("Working in Life Sciences") xlabel(-3 "-7" -2 "-5" -1 "-3"  0 "-1"  1 "1" 2 "3" 3 "5", labsize(small))  legend(order(1 "TWFE OLS" 3 "Callaway-Sant'Anna") rows(1) region(style(none))) xline(0, lcolor(gs8%60) lpattern(dash)) yline(0, lcolor(gs8%60) lpattern(dash)) graphregion(color(white)) bgcolor(white) ylabel(, angle(horizontal))) lag_opt1(msymbol(Oh) color(dknavy)  msize(small)) lag_ci_opt1(color(dknavy%100)) lag_opt2(msymbol(Dh) color(cranberry) msize(small)) lag_ci_opt2(color(cranberry)) lag_opt5(msymbol(Sh)  msize(large)  color(black)) lag_ci_opt5(color(black)) perturb(-0.05(0.07)0.01) alpha(0.05)
graph export "${path}\output\figures\acs-more_twfecs_nevert.pdf", replace